#include "FracSolver.hpp"

#ifdef _WIN64
const std::string data_path_prefix = "D:/HPCA/OpenAtom_SolAsymDense/StepMatrixRHSAndSolution-3";
#else
const std::string data_path_prefix = "./StepMatrixRHSAndSolution-1";
#endif

int main(int argc, char* argv[])
{
    int arg_cnt = 1;
    int solver_type = atoi(argv[arg_cnt++]);
    int prec_type = atoi(argv[arg_cnt++]);
    int firstId = atoi(argv[arg_cnt++]);
    int finalId = atoi(argv[arg_cnt++]);
    // SequenceSolver solver(solver_type, prec_type, firstId);
    FracSolver solver(solver_type, prec_type, 128);

    int test_cnt = 5;
    printf("Test %d times\n", test_cnt);

    for (int case_id = firstId; case_id <= finalId; case_id++) {
        char buf[100];
#ifdef _WIN64
        sprintf_s(buf, "%s/%05d", data_path_prefix.c_str(), case_id);
#else
        sprintf(buf, "%s/%05d", data_path_prefix.c_str(), case_id);
#endif
        std::string file_prefix(buf);

        int nrows, ncols;
        double* mat = nullptr, * rhs = nullptr, * exact_sol = nullptr;
        read_Fortran_unformatted((file_prefix + ".rhs"), nrows, ncols, rhs);        // check_sum(nrows, rhs, 10);
        read_Fortran_unformatted((file_prefix + ".sol"), nrows, ncols, exact_sol);  // check_sum(nrows, exact_sol, 10);
        read_Fortran_unformatted((file_prefix + ".matr"), nrows, ncols, mat);       // check_sum(nrows * ncols, mat, 10);
        assert(nrows == ncols);
        double* tmp_mat = new double [nrows * nrows];
        double* tmp_rhs = new double [nrows];
        // 开辟一些没什么用的数组，全部填0，用以5个分块的处理
        std::vector<double> dummy_A13( (nrows / 2) * (nrows / 2), 0.0 ),
                            dummy_A23( (nrows / 2) * (nrows / 2), 0.0 ),
                            dummy_A33( (nrows / 2) * (nrows / 2), 0.0 ),
                            dummy_A32( (nrows / 2) * (nrows / 2), 0.0 ),
                            dummy_A31( (nrows / 2) * (nrows / 2), 0.0 );
        std::vector<double> dummy_x3 ( nrows / 2, 0.0 ),
                            sol      ( nrows    , 0.0 ), 
                            dummy_b3 ( nrows / 2, 0.0 );

        for (int t = 0; t < test_cnt; t ++) {
            copy_vec(nrows * nrows, mat, tmp_mat);
            copy_vec(nrows        , rhs, tmp_rhs);
            
            solver.solve(nrows / 2, nrows / 2, tmp_mat,
                dummy_A13.data(), dummy_A23.data(), dummy_A31.data(), dummy_A32.data(), dummy_A33.data(),
                tmp_rhs, tmp_rhs + nrows / 2, dummy_x3.data(), sol.data(), dummy_b3.data());
        }

        if (check_correctness(nrows, mat, sol.data(), rhs)) printf("Accept!\n");
        else printf("Reject.\n");

        write_Fortran_unformatted(file_prefix + ".submitSol", nrows, 1, sol.data());

        delete[] mat; delete[] rhs; delete[] exact_sol;
        delete[] tmp_mat; delete[] tmp_rhs;
    }

    return 0;
}
